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Abstract 

We derive explicit formulas for the reconstruction of a function from its integrals 
over a family of spheres, or for the inversion of the spherical mean Radon trans- 
form. Such formulas are important for problems of thermo- and photo- acoustic 
tomography. A closed-form inversion formula of a filtration-backprojection type is 
found for the case when the centers of the integration spheres lie on a sphere in 

surrounding the support of the unknown function. An explicit series solution is 
presented for the case when the centers of the integration spheres lie on a general 
closed surface. 



Introduction 

The problem of the reconstruction of a function from its spherical integrals (or means) 
have recently attracted attention of researchers due to its connection to the thermo- 
acoustic and photo-acoustic tomography [12,13,21,22]. In these imaging modalities, the 
object of interest is illuminated by a short electromagnetic pulse which causes a fast 
expansion of the tissue. The intensity of the resulting ultrasound wave is recorded by 
a set of detectors surrounding the object. The local intensity of such expansion is of 
significant medical interest: it depends on the physical properties of the tissue (such as, 
for example, water content) and its anomaly can be indicative of tumors. Under certain 
simplifying assumptions the measurements can be represented by the integrals of the 
expansion intensity over the spheres with the centers at the detectors locations. The 
reconstruction of the local properties from these integrals is equivalent to the inversion of 
the spherical mean Radon transform. 

An introduction to the subject can be found in [12-14]; for the important results on 
the injectivity of the spherical Radon transform and the corresponding range conditions 
we refer the reader to [2-4, 9, 10]. In the present paper we concentrate on explicit inversion 
formulas which are important from both theoretical and practical points of view. Most of 
the known formulas of this sort pertain to the spherical acquisition geometry, i.e. to the 
situation when centers of the integration spheres (the positions of the detectors) lie on 
a sphere surrounding the body. Such are the series solutions for 2-D and 3-D presented 
in [15,16,21]. More desirable backprojection-type formulas were derived in [9] for odd- 
dimensional spaces, and implemented in [6] . A different explicit formula for the spherical 
acquisition geometry valid in 3-D was found in [22] (together with formulas for certain 
unbounded acquisition surfaces). 

In this paper we present a set of closed-form inversion formulas for the spherical ge- 
ometry and a series solution for certain other measuring surfaces. Our formulas for the 
spherical case are of the filtration-backprojection type; they are valid in M", n > 2. 
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Such formulas for the even- dimensional cases were not known previously. (A set of dif- 
ferent inversion formulas for the even- dimensional case was announced by D. Finch [11] 
during tomography meeting in Oberwolfach in August, 2006, where our results where also 
presented for the first time.) 

The spherical case is discussed in Sections 1 through 3. A series solution for a general 
acquisition geometry is presented in Section 4. 

1 Formulation of the problem 

Suppose that Cq function /(x), x gM", n > 2 is compactly supported within the closed 
ball B of radius R centered at the origin. We will denote the boundary of the ball by 
dB. Our goal is to reconstruct /(x) from its projections g{z, r) defined as the integrals of 
/(x) over the spheres of radius r centered at z: 

t,(z,r)= j fiz + riy-^dsii). 

gn-l 

where §"~^ is the unit sphere in M", £ is a unit vector, and ds is the normalized measure 
in M". Projections are assumed to be known for all z G dB, < r < 2R (integrals 
for r > 2R automatically equal zero, since the corresponding integration spheres do not 
intersect the support of the function). In the following two sections we will present an 
explicit formula of backprojection type that solves this reconstruction problem. 



2 Derivation 

Our derivation is based on certain properties of the solutions of the Helmholtz equation 



m 



A/i(x) + A2/i(x) = 0. 



For this equation the free space Green's function $(x, y. A) satisfying radiation boundary 
condition is described (see, for example [1]) by the formula 



. / N n/2-1 

where H^^^_^{t) is the Hankel function of the first kind and of order n/2 — 1. To simphfy the 
notation, we introduce functions J{t), N{t), and H{t), defined by the following formulas 



J{t) 

m 



Jn/2-l{'t) 
^n/2-1 ' 

Nn/2-l{t) 

r(i) 



H{t)- -Jit) + ^N{t), 
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where Jn/2-i{t) and A^„/2-i(^) are respectively the Bessel and Neumann functions of order 
n/2 — 1. In this notation Green's function $(x, y, A) can be re-written in a simpler form: 



$(x,y,A) = ic(A,n)//(A|x-y|) 

= c(A,n)[iJ(A|x-y|)-Ar(A|x-y| 

where c(A, n) is a constant for a fixed value of A: 

\n-2 

c(A, n) 



4(2vr: 



|n/2-l ' 



We note that function J(A|x|) is a solution of the Helmholtz equation for all x e M", 
while A^(A|x|) solves this equation in R"\{0}. 

In other to derive the inversion formula, we utilize the following integral representation 
of / in the form of a convolution with J(A|y — x|): 



/(y) = ^^^ / (y"/(x)J(A|y-x|)cia;j A-^dA. 

The above equation easily follows from the Fourier representation of / 
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(1) 



(27r)" 
1 

(27r)" 




e-'--V(x)dxde 




/(x) 



and from the well-known integral representation for J(|u|) [17] 

gn-l 

Let us denote the inner integral in (1) by G'j(y, A) 



Gj{y,X)= J /(x)J(A|y-x|)cix. 



(2) 



Similarly to the kernel J(A|y — x|) of this convolution, function G'j(y, A) is an entire 
solution of the Helmholtz equation. Boundary values of this function Gj{z,X), z G dB 
are easily computable from projections: 



2R 



Gj{z,X)= / /(x)Jo(A|z-x|)(ix= / Jo{Xr)g{z,r)dr. 



If A is not in the spectrum of the Dirichlet Laplacian on B, Gj(y, A) is completely deter- 
mined by its boundary values and can be found by solving numerically the corresponding 
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Dirichlet problem for the Helmholtz equation. Then /(x) can be reconstructed from 
equation (1). Unfortunately, such a solution would not have an explicit form. 

In order to obtain an explicit formula for Gj(y, A) we will utilize a Helmholtz repre- 
sentation for J(A|y — x|); it results from an application of Green's formula and has the 
following form: 



J(A|y-x|) = J 

dB 



J(A|z- 



d 



I — 
dn. 



■*(y,z,A) 



-$(y,z,A)^J(A|z-x| 



or 



J(A|y -x|) = -c(A,n 

dB 

-N{X\y-z\ 



J(A|z-x 

d 



_d_ 
dn. 



J(A|z-x|) 



iV(A|y-z|) 



ds{z). 



(3) 



Such a representation is valid for any bounded single-connected domain with sufficiently 
regular boundary. A straightforward substitution of equation (3) into (2) leads to a 
boundary value representation for Gj{y, A) involving the normal derivative the latter 
function. Unlike the boundary values of G'j(z, A), the normal derivative ^Gj(z, A) 
cannot be explicitly computed from projections g{z,r). 

This difficulty can be circumvented by modifying the Helmholtz representation as 
described below. We notice that in the special case of a spherical domain the second 
integral in the formula (3) 



7(x,y) = J 7V(A|y-z|)Aj(A|z-x|Ms(z) 



(4) 



dB 



is a symmetric function of its arguments i.e. that 

/(x,y) = /(y,x). 

The proof of this fact is presented in the Appendix. Using this symmetry we obtain a 
modified Helmholtz representation for J(A|y — x|): 



J(A|y-x|) = -c(A,n) J 

dB 

-Ar(A|z-x|) 



J(A|z-x| 



_d_ 



.Ar(A|y-z|) 



ds{z). 



(5) 



Now the substitution of equation (5) into (2) yields 
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J /(x) J(A|y - ^\)dx = -c(A, n)J I J /(x) J(A|z - ^\)dx j -^^N{X\y - z|) 

B dB \.\b / 

J /(x)7V(A|z-x|)c^xJ ^^(A|y-z|) 



ds{z), 



where the inner integrals are easily computable from the projections g{z,t): 



2R 



J /(x)J(A|z-x|)cix = J J{Xt)g{z,t)dt, 





2R 



f{x)N{X\z-x\)dx = j N{Xt)g{z,t)dt. 

B 

Thus, the convolution of / and J can be reconstructed from the projections as follows: 
J /(x) J(A|y - ^\)dx = -c(A, n) J U J{Xt)g{z, t)dt \ ^^N{X\y - z|) 



dB 



r 2R \ 
j N{Xt)g{z,t)dt\ Aj(A|y-z|) 



ds{z) 



c(A, n)div j n(z) 

dB 



2R 



J{Xt)g{z,t)dt\N{X\y-z\) 



2R 



N{Xt)g{z,t)dt\ J(A|y-z| 



ds{z). 



(6) 



Finally, by combining equations (1) and (6), one arrives at the following inversion formula 

/(y) = 4(27rV-^ ^^^ / '^(^)^(^'ly 

where 



dB 



h{z,t) 



2R 



N{Xt) I J J{Xt')g{z,t')dt' 



2R 



J(Xt) J N(Xt')g{z,t')dt' 



X'^-^dX. 



3 Particular cases 



3.1 2-D case 



From the point of view of practical applications the two- and three- dimensional cases are 
the most important ones. In 2-D, J{t) = Jo{t), N{t) = NQ(t), and the inversion formula 
has the following form 



/(y) = ^div j n(z)/i(z, |y - z\)dl(z), 

dB 



(7) 



where 



h{z,t) 



2R 

NoiXt) I I MXt')giz,t')dt' 

2R 

Jo{Xt) I J No{Xt')g{z,t')dt' 



XdX. 



(8) 



Equation (7) is a backprojection followed by the divergence operator. It is worth noticing 
that such a divergence form of a reconstruction formula is not unusual; it also naturally 
occurs in reconstruction formulas for the attenuated Radon transform (see, for instance 
[20]). 

Equation (8) represents the filtration step of the algorithm. In order to better under- 
stand the nature of this operator we re-write (8) in the form 



h{z,t) 



2R 



Hi^\Xt) { j H^^\Xt')g{z,t')dt' 



2R 



-Hl,'\Xt)lj4'\Xt')g{z,t')dt' 



XdX. 



(9) 



and recall that for large values of the argument t the Hankel function HQ^\t) has the 
following asymptotic expansion ([19]) 

In the situation when the support of the function /(x) remains bounded and the radius R 
of the ball B becomes large, the inner and outer integrals in (9) reduce to the direct and 
inverse Fourier transforms, with one of the terms corresponding to the positive frequencies 
and the other to the negative ones. Thus, since the two terms have opposite signs, in the 
asymptotic limit of large R operator (9) equals (up to a constant factor) to the Hilbert 
transform of 5f(z, •). 
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3.2 3-D case 

In the three-dimensional case 



Xr 

Ni/2{Xr) 



Xr 



Jo(Ar), 



-no(Ar), 



where jo{t) and no{t) are the spherical Bessel and Neumann functions respectively. The 
formula takes the following form 



/(y)=Y^div y n(z)/i(z, |y-z|)cis(z), 



(10) 



\z\=R 



with 



2R 



no{Xt) / jo{Xt')g{z,t')dt' 



2R 



-jo{>^t) U no{Xt')g{z,t')dt' 



X^dX. 



In this case, however, a further simplification is possible since jo{t) and no{t) have a simple 
representation in terms of trigonometric functions: 

. , . sini , , cost 
Jo{t) = — — , no{t) = 



t 



t 



The substitution of these trigonometric expressions into the inversion formula leads to a 
significantly simpler formula: 



2 f 

h{z,t) = / cos(At) 

TTt J 



M+ 



2R 



XdX 



+ — [ sin(At) 
nt J 



2R 



cos{Xt')^^^dt' 



t' 



XdX 



2R 



TTt dt 



cos(At) 

'2d g{z,t) 
tdt t ' 



.0 

2R 



dX 



J cos(At) J cos{Xt')^^^dt' 



dX 



(11) 
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where we took into account the fact that the Fourier sine and cosine transforms are self- 
invertible. By combining (11) and (10) our inversion formula can be re-written in the 
form 



This expression is equivalent to one of the formulas derived in [22] for the 3-D case. 

4 Inversion of the spherical mean Radon transform 
in other geometries 

In such applications as photo- and thermo- acoustic tomography, the designer of the 
measuring system has a freedom in selecting the detectors' locations. The detectors (the 
centers of the integration spheres) do not have to lie on a sphere. In this section we 
present reconstruction formulas for the case when the measuring surface is a boundary of 
certain other domains. Namely, our method works for the domains whose eigenfunctions 
Urai'x) of the (zero) Dirichlet Laplacian are explicitly known. Such are, for example, the 
domains for which the eigenfunctions can be found by the separation of variables, i.e. 
sphere, annulus, cube, and certain subsets of those, and crystallographic domains (see 



The proof of the range theorem in [2] involves implicitly a reconstruction procedure 
also based on eigenfunction expansions. Unlike the present method, that procedure would 
involve division of analytic functions that have countable number of zeros. While the range 
theorem guarantees cancellation of these zeros when the data are in the range of the direct 
transform, a stable numerical implementation of such division would be complicated if 
not impossible. (Similarly, the scries solution of [15] for 2-D circular geometry involves 
division by Bessel functions.) The technique we present below does not require such 
divisions. 

Suppose A^, Um{x.) are the eigenvalues and eigenfunctions of the (negative) Dirichlet 
Laplacian on a bounded domain Q with zero boundary conditions, i.e. 



As in the previous sections, we would like to reconstruct a function /(x) G L^(Q,) from 
the known values of its spherical integrals g{z,r) with the centers on dfl: 



We notice that Wm(x) is a solution of the Dirichlet problem for the Helmholtz equation 




[7,8]). 



Am„(x) + A>^(x) = 0, ^en, new, 

Mto(x) = 0, x e dfl. 




z e dQ. 
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and thus admits the Helmholtz representation 

f 9 
Jan '^^ 



f d 

^ic{Xm,n) H{\m\x- z\)—Um{z)ds{z) xeQ. (12) 

Jan 

On the other hand, eigenfunctions {umi'x.)}'^ form an orthonormal basis in L2(il). 
Therefore /(x) can be represented (in L^sense) by the series 

oo 

/(x) = 2^ amUmi^) (13) 
m=0 

with 

otm= I Umi?^)f{'^)d:i^. (14) 

The reconstruction formula will result if we substitute representation (12) into (14) and 
change the order of integrations 

ttm = / 'Um(x)/(x)(ix 

^ic{Xm,n) [ ([ if(A„|x-z|)^M„(z)ds(z) ) /(x)dx 
Jn \Jdn / 

^ic{Xm,n)j (^j H{Xm\x-z\)f{x)dxj -^Um{z)ds{z). (15) 

The change of the integration order is justified by the fact that eigenfunctions Mm(x) are 
continuous. The inner integral in (15) is easily computed from projections 

/ if(A„,|x-z|)/(x)dx= / g{z,r)H{Xmr)dr, 
Ja J 

SO that 

With Fourier coefficients now known /(x)is reconstructed by summing series (13). 
If desired, this solution can be re- written in the form of a backprojection-type formula: 

/(x) = ^ q;^w^(x) = / i^amc{Xm,n)iH{Xm\y^-z\)—u„^{z)\ds{z) 

m=0 "^^^ \m=0 ^ / 

= / /i(z, |x-z|)ds(z), (17) 
Jan 

where 

h{z, t)^i^ c{Xm, n)amH{Xmt)-^Um{z), (18) 



dn 

m=0 
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and coefficients are computed using equation (16). In the above formula equation (17) 
is clearly a backprojection operator, and (18) is a filtration. However, the latter operator 
is now represented by a series rather than by a closed form expression. 

The series solution described above have an interesting property not possessed (to the 
best of our knowledge) by any other explicit reconstruction technique. Let us consider a 
slightly more general problem. Suppose that region f2 is a proper subset of a larger region 
fli {Q C Qi) and that a function F is defined on fli. We will denote the restriction of 
F on O by /, i.e. 

' F(x) , X e Q 
, X e w\n ■ 



/(x) 



We would like to reconstruct /(x) from the integrals g{z,r) of F over spheres with the 
centers on dQ: 

^(z,r)= J F(z + rs)r"-^ds, z e dQ. 

§n-l 

The difference with the previously considered problem is in that the centers of the inte- 
gration spheres are know lying on a surface which is inside the support Qi of the function 
F. While we are still trying to reconstruct the restriction / of F to ^2, the integrals we 
know are those of F and not of /. 

It turns out that the solution to this problem is still given by formulas (13) and (16) 
(or equivalently by (17), (18), and (16)). Indeed, if we extend functions Mm(x) by to 
M"\r2, formula (12) holds for all x G M**, and (13) remains unchanged. In formula (15) / 
can be replaced by F as follows: 



M^(x)/(x)rfx = / M„(x)F(x)(ix 
Jn 

^ic{Xm,n) j {^j if(A^|x- z|)F(x)(ix^ ^M„(z)ds(z). 



and the inner integral can be computed from projections as before: 

i/(A„|x-z|)F(x)cix= / g{7^,r)H{\mT)dr. 



/ 



By combining the two above equations we again arrive at the formula (16). 

To summarize, if the centers of the integration spheres lie on a closed surface inside 
of the support of the function, the values of the function corresponding to the interior of 
that surface can be stably reconstructed by a formula containing the normal derivatives 
of the eigenfunctions of the Dirichlet Laplacian. If these eigenfunctions are given by an 
explicit formula, our technique yields an explicit series solution to this problem. 



Appendix 

In this section we prove that for arbitrary x, y G i? function /(x, y) defined by equation (4) 
is a symmetric function of its arguments, i.e. that /(x, y) = /(y,x). 
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First we need to find a closed form expression for a single layer potential of a spherical 
shell with the density equal to a spherical harmonic Y^^x). To this end we introduce 
functions m^''(x) and ^^''(x) defined as follows 

x.5''(x)=l^(^)(x)J(,)(A|x|), 
where we, as before, utilize the notation 

Jn/2+k-l{t) 



tn/2-1 



n/2+k-l 



it) 



H{k){t) 



^n/2-l • 
^n/2-1 ' 



Functions ti^''(x) are ([18]) entire solutions of the Helmholtz equation for all x e IR'*, 
while v'l'\'x.) are radiating solutions of this equation in R"\{0}. If we apply the Green's 
theorem in an exterior of a sphere of radius ro to functions (x) and v^' (x), we will 
obtain (for |x| > ro) 



«a''(z)^*(x, z. A) - $(x, z, X)^u'^\z) 



dz 



|z|=ro 



Y,^'\z) ^J(fe)(Aro)^$(x,z, A) - $(x, z, A)AJ(',)(Aro) 



dz, 



\z\=ro 



and 



k,l / \ 

(x) 



^A''(^)^^(x, z, A) - $(x, z. A) A^M(z) 



dz 



|z|=ro 



d 



^^■(fe)(Aro)^$(x,z, A) - $(x,z, A)Aii'('fc)(Aro) 



dz 



|z|=ro 



where 5 = z/|z|. By combining the above two equations one can eliminate the terms with 
a . 
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J{k){\rQ)v^ (x) = A [i/(fc)(Aro)J(fc)(Aro) - J(fc)(Aro)i/(fc)(Aro)] 
X j Y'/''^(5)$(x,z,A)dz. 



(19) 



|z|=ro 



The expression in the brackets can be simplified using the well known formula ([19]) for 
the Wronskian of Ha\t) and Ja{t)- 
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27ri 



Formula (19) then takes form 

|z|=ro 

SO that the single layer potential we consider is described by the following equation 

I r/'H;S)$(x,z, A)dz = z27rro(Aro)'^-V(fe)(Aro)%)(A|x|)F/')(£). (20) 

|z|=ro 

By substituting into (20) the expression for the Green's function $(x, z, A) in the form 

$(x, z, A) = ic(A, n)i^(A|x - z|) 

one obtains 

I yI'\z)H{\\^ - z|)rfz = '^'''f^''^^'' Jik) (Aro)%) (A|x| )y^(^) (x) . 

\z\=ro 

Finally, by separating the real and imaginary parts of the above equation we arrive at the 
following two formulas: 

I F/'^\z)J(A|x-z|)dz = ^^^^^^^ (21) 



jz|=ro 



I Y,^'\z)NiM^ - ^\)dz = ^""'.(A^^f" J(,)(Aro)7V(,)(A|x|)l^(^^(x), (22) 

M=ro 

valid for |x| > tq. This completes the preparation for the proof of the symmetry of /(x, y). 
Let us consider function I{ax,(3y). For fixed values of a and (3 this is an infinitely 

smooth function of x and y defined on S"~^ x S"~^. Consider the Fourier expansion of 
I{ax, I3y) in the spherical harmonics in both variables x and y. Coefficients of such series 
are given by the formula 

/ / Y,\x)Yp\y)I{ax,(3y)dxdy (23) 



an-l Sn-1 



with fc = 0,1,2,..., 0</<4, 

fn + k-l\ fn + k-3\ , ^ 

di = n,do = 1. 

In order to prove the symmetry 7(x, y) = 7(y, x) it is enough to prove that 

= %f (A") (24) 
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for all relevant values of k, k', I, I'. By substituting the expression for I{ax, (5y) into (23) 
one obtains 



%n-l sn-1 



d 

iV(A|z - ax\)^J{(5y - z|)dz 



|z|=R 



dxdy 



yv\y) 



a 



n—l 



k'l 



Yi\x)N{X\z - ax\)dx 



_d_ 



J{X\f3y-z\dz 



Yi\u/\u\)N{X\z-u\)du 



Vl=|x| 



/ 



_d_ 



J{X\(3y-z\dz 



Utilizing formulas (21), (22) and the orthonormality of the spherical harmonics we find 



that al',li,{a, P) = ii k ^ k' or I ^ I'. Otherwise 



c(A, n) 



i?^- ( Aa) J(fc) ( A/5)iV(fe) ( Ai?) J(',) ( Ai?) . 



Inspection of the above formula shows that coefficients a^',li,{a, (5) indeed satisfy (24) and, 
thus, that /(x, y) = 7(y, x). 
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